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Abstract 

We present a new cosmoiogicai gaiaxy formation modei, v^GC, as an updated version of our 
previous modei i^GC. We adopt the so-caiied “semi-anaiytic” approach, in which the formation 
history of dark matter haios is computed by A^-body simuiations, whiie the baryon physics such 
as gas cooiing, star formation and supernova feedback are simpiy modeied by phenomenoiog- 
icai equations. Major updates of the modei are as foiiows; (1) the merger trees of dark matter 
haios are constructed in state-of-the-art A^-body simuiations, (2) we introduce the formation 
and evoiution process of supermassive biack hoies and the suppression of gas cooiing due 
to active gaiactic nucieus (AGN) activity, (3) we inciude heating of the intergaiactic gas by the 
cosmic UV background, and (4) we tune some free parameters reiated to the astrophysicai 
processes using a Markov chain Monte Cario method. Our A^-body simuiations of dark matter 
haios have unprecedented box size and mass resoiution (the iargest simuiation contains 550 
biiiion particies in a 1.12 Gpc/h box), enabiing the study of much smaiier and rarer objects. 
The modei was tuned to fit the iuminosity functions of iocai gaiaxies and mass function of neu- 
trai hydrogen. Locai observations, such as the Tuiiy-Fisher reiation, size-magnitude reiation 
of spirai gaiaxies and scaiing reiation between the buige mass and biack hoie mass were weii 
reproduced by the modei. Moreover, the modei aiso weii reproduced the cosmic star formation 
history and the redshift evoiution of rest-frame TT-band iuminosity functions. The numericai 
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catalog of the simulated galaxies and AGNs is publicly available on the web. 
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1 Introduction 

Understanding the formation and evolution of galaxies is a pri¬ 
mary goal in astrophysics. Over the past decades, wide and deep 
surveys at various wavelengths have acquired numerous obser¬ 
vational data of galaxies spanning a wide range of galaxy types, 
magnitudes and distances (see Madau & Dickinson 2014, for re¬ 
view). Theoretically, the A cold dark matter (CDM) paradigm 
can explain the formation of the large scale structures governed 
by dark matter (DM) and dark energy. However, at the scale of 
galaxies, where baryons play important roles, several inconsis¬ 
tencies remain between the theory and observations. To fully 
elucidate galaxy formation, we need to solve the complicated 
physical processes of baryons within the framework of A-CDM 
universe. 

One of the most promising ways to address this issue is the 
hydrodynamical simulations of cosmological galaxy formation, 
in which the equations of gravity, hydrodynamics, and thermo¬ 
dynamics are solved self-consistently. However, the mass reso¬ 
lution and box size of these simulations are still limited by com¬ 
putational costs, and the physical processes on scales smaller 
than the numerical resolution are treated by phenomenological 
recipes (the so-called “sub-grid physics”), which contain large 
uncertainties (see Springel 2012, for review). 

“Semi-analytic models” (SA models) are also widely used 
in studies of cosmological galaxy formation (e.g., Kauffmann 
et al. 1993b; Cole et al. 1994, 2000; Somerville & Primack 
1999). In SA models, the formation and evolution history of 
dark matter halos are explicitly modeled by analytical formulae 
or A-body simulations, while the complicated baryon physics 
are modeled by phenomenological equations. The advantage 
of this technique is its lower computational cost than numer¬ 
ical simulation, enabling us to create a large sample of mock 
galaxies covering the wide range of physical properties such as 
mass, magnitude, and spatial scale. We can also investigate a 
wide range of the parameter space and test various models of 
the baryon physics. However, to discuss the galaxy-scale dy¬ 
namics, we need to combine SA models (which do not explic¬ 
itly treat such dynamics) with fully numerical simulations. See, 
e.g., Somerville & Dave (2014) for more detailed review of the 
physical models of cosmological galaxy formation. 

In this paper we introduce our new galaxy formation model. 
New Numerical Galaxy Catalog (r/^GC), an updated version of 
Numerical Galaxy Catalog {vGC) presented in Nagashima et al. 
(2005, hereafter N05; see also Nagashima & Yoshii 2004). Our 
model is an SA model, in which we directly extract the merger 
trees of DM halos from A-body simulations, following the pio¬ 


neering work of Roukema et al. (1997). The vGC model and its 
variants have been used in many studies (e.g., Kobayashi et al. 
2007, 2010; Okoshi et al. 2010; Makiya et al. 2011, 2014; Enoki 
et al. 2014; Shirakata et al. 2015; Oogi et al. 2016). Major up¬ 
dates of the new model from the version of N05 are as follows: 
(1) u^GC adopts the new A-body simulations of DM halos re¬ 
cently presented by Ishiyama et al. (2015), (2) the formation 
and evolution process of supermassive black holes (SMBHs) 
and suppression of gas cooling by active galactic nuclei (AGNs) 
are included, (3) heating of the intergalactic gas by the cosmic 
UV background is included, and (4) some parameters are tuned 
to lit the local luminosity functions and HI mass function using 
a Markov chain Monte Carlo (MCMC) method. 

Several other groups have also proposed SA models (see 
Somerville & Dave 2014, for review). Each of these models 
is based on different A-body simulations and adopts different 
equations of the baryon physics. For a comparison study of dif¬ 
ferent galaxy formation models, see Knebe et al. (2015). Our 
model is characterized by the substantially higher mass resolu¬ 
tion of the A-body simulations of DM halos, comparing with 
other large box simulations. Our simulations consist of seven 
runs with varying mass resolutions and box sizes, as listed in 
table 1. For example, the largest simulation, v^GC-L, includes 
8192® DM particles in a box of 1.12 h~^ Gpc, and the mini¬ 
mum halo mass reaches 8.79 x 10 ®Mq. Comparing with the 
Millennium simulation (Springel et al. 2005), the i/^GC-L sim¬ 
ulation is four times better in mass resolution and is 11 times 
larger in spacial volume. The o^GC-Hl simulation has the 
highest mass resolution among our simulations. The minimum 
halo mass reaches 1.37 x 10® Mq, below the effective Jeans 
mass at high redshift (N05). This mass resolution is 2 times bet¬ 
ter than Millennium-II simulation (Boylan-Kolchin et al. 2009), 
although the spatial volume of o^GC-Hl is 3 times smaller than 
that of Millennium-II. These high mass resolution and large spa¬ 
tial volume enable us to obtain a statistically significant number 
of mock galaxies and AGNs, even at high redshifts. Moreover, 
we adopt the cosmological parameters recently obtained by the 
Planck satellite (Planck Collaboration et al. 2014), while the 
most of other SA models are based on the parameters obtained 
by Wilkinson microwave anisotropy probe (WMAP), which sig¬ 
nificantly differ from the Planck results. For more detailed 
comparison with other cosmological A-body simulations, see 
Ishiyama et al. (2015). 

This paper describes the basic properties of our model, fo¬ 
cusing on the nature of local galaxies. The properties of distant 
galaxies and AGNs will be discussed in our forthcoming papers. 
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The paper is organized as follows. Sections 2 and 3 present 
the details of our model and the parameter fitting method, re¬ 
spectively. The general properties of our numerical galaxy cat¬ 
alog are presented in section 4, and sections 5 and 6 compare the 
model predictions with the observed properties of local and dis¬ 
tant galaxies. Section 7 summarizes the paper. The mock galaxy 
catalog produced by the our new model is publicly available on 
the web'. 


2 Model descriptions 

In the CDM universe, DM halos hierarchically grow from small 
to large scales. When a DM halo collapses, the contained gas is 
heated to virial temperature by shock, and then gradually cools 
by radiative cooling (in reality, a gas in low-mass halos would 
not he shock heated but directly forms cold gas disk; see sec¬ 
tion 2.2 for more detailed discussion). The cooled gas condense 
into stars; those stars and dense cold gas constitute galaxies. 
The massive stars formed hy this process explode as supernovae 
(SNe), blowing out surrounding cold gas. This process sup¬ 
presses further star formation (the so-called “SN feedback”). 
Massive stars also eject metals. Galaxies in a common DM halo 
sometimes merge into more massive galaxies, and galaxy bulge 
is formed as a merger remnant; cold gas in the merger rem¬ 
nant is converted into stars with short timescale, a phenomenon 
called a starburst. During the starburst, a fraction of the cold 
gas gets accreted by the supermassive black hole (SMBH) at 
the galaxy center. By repeating these processes, galaxies and 
SMBHs have formed and evolved to the present epoch. Each 
of these processes is described in the following subsections. 
Figure 1 displays an overview of the model. 

2.1 Dark matter merger trees 

The merger trees of DM halos are directly extracted from a se¬ 
ries of large cosmological Mbody simulations, called the v^GC 
simulations (Ishiyama et al. 2015). The basic properties of the 
v^GC simulations are summarized in table 1. We conducted 
seven simulations, varying the mass resolution and spatial vol¬ 
ume. The largest u^GC-L run simulated the motions of 8192® 
(550 billion) DM particles in a comoving box of 1.12 h“'Gpc. 
The mass resolution was 2.20 x 10® which is the best 

one among simulations applying boxes larger than l/i“'Gpc. 
The mass resolution of the run with the smallest box (i^^GC- 
H2) was 3.44 x 10® which is sufficient to resolve small 

dwarf galaxies. By combining these simulations, we can gener¬ 
ate mock catalogs of galaxies and AGNs with unprecedentedly 
high resolution and statistical power. 

The cosmological parameters of the u^GC simulations were 
based on the concordance ACDM model consistent with ob- 

‘ http://www.imit.chiba-u.jp/faculty/nngc/ 


servational results obtained by the Planck satellite (Planck 
Collaboration et al. 2014). Namely, Do = 0.31, Db = 0.048, 
Ao = 0.69, h — 0.68, Us = 0.96, and as = 0.83. The u^GC sim¬ 
ulations were conducted by using a massively parallel TreePM 
code GreeM (Ishiyama et al. 2009, 2012). DM halos are iden¬ 
tified by the friends-of-friends (FoF) group finder (Davis et al. 
1985), with the linking parameter b = 0.2. The smallest ha¬ 
los consisted of 40 particles. The spatial positions of the halos 
were tracked by using those of the most bound particles. It has 
been known that the properties of halo merger tree are depend 
on the halo finding algorithm and tree building algorithm (see, 
e.g., Knebe et al. 2011; Onions et al. 2012; Elahi et al. 2013; 
Srisawat et al. 2013; Avila et al. 2014; Lee et al. 2014). For 
further details of the u^GC simulations and the method for ex¬ 
tracting the merger trees, see companion paper, Ishiyama et al. 
(2015). 


2.2 Gas cooling 

We define the formation epoch of the DM halo as the time at 
which the DM halo mass doubles its mass since the last for¬ 
mation epoch (Lacey & Cole 1993). At this time, the physical 
quantities of the halo, such as circular velocity, halo age and 
mass density are re-estimated. Before reionization of the uni¬ 
verse, the mass fraction of the baryonic matter in a collapsing 
DM halo is given by (/b) = Db/Dp (after cosmic reionization, 
the baryon mass in a halo deviates from (/b); see subsection 
2.3). The baryonic matter consists of diffuse hot gas, dense cold 
gas, stars, and black holes. When a mass of DM halo decreases, 
diffuse hot gas also decreases at the same ratio with the decrease 
of DM mass, while a mass of other baryon components does not 
change. 

When a DM halo of circular velocity Kirc forms, the con¬ 
tained gas is shock heated to the virial temperature Tvir of the 
halo: 


Tvi 


1 pnip 

2 ks 


1/2 
circ ? 


( 1 ) 


where mp,k^ and p are the proton mass, Boltzmann constant 
and mean molecular weight, respectively. Following Shimizu 
et al. (2002), we assume that the hot gas distributes through the 
DM halo with an isothermal density profile wifh a finite core 
radius: 


Phot{r) 


PhotjO 

1 + (r/rc^ ’ 


( 2 ) 


where Vc = 0.227?vir/c, and Rvu is the virial radius of the host 
halo. The concentration parameter c is known to be a function of 
DM halo mass and redshift. We used the analytical formula of 
c proposed by Prada et al. (2012), which is obtained by fitting 
cosmological A/-body simulations. The model of Prada et al. 
(2012) and Sanchez-Conde & Prada 2014 are consistent with 
our M-body simulations. 
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Fig.l. Schematics of the model. (Left) Flow chart of the model showing how the model predicts the observable properties of galaxies. (Right) Schematic of 
the transfer of baryon components. 

Table 1. Details of the u^GC simulations. N is the number of simulated particles, L is the comoving box size, m is the particle mass 
resolution, Mmin is the mass of the smallest halos, the total number of halos, and Mmax is the mass of the largest halo in each 
simulations. The smallest halos consist of 40 particles. In the last two columns, values at z=0 are presented except for the 
simulation, which was stopped at ^ = 4._ 


Name 

N 

L(h-®Mpc) 

m{h ^Mq) 


#Halos 

Afinax(h ^Mq) 

i/^GC-L 

8192® 

1120.0 

2.20 X 10® 

8.79 X 10® 

421,801,565 

4.11 X 10®® 

i^^GC-M 

4096® 

560.0 

2.20 X 10® 

8.79 X 10® 

52,701,925 

2.67 X 10®® 

u^GC-S 

2048® 

280.0 

2.20 X 10® 

8.79 X 10® 

6,575,486 

1.56 X 10®® 

ix^GC-SS 

512® 

70.0 

2.20 X 10® 

8.79 X 10® 

103,630 

6.58 X 10®"® 

jx^GC-HI 

2048® 

140.0 

2.75 X 10^ 

1.10 X 10® 

5,467,200 

4.81 X 10®'® 

u'^GC-m 

2048® 

70.0 

3.44 X 10® 

1.37 X 10® 

4,600,746 

4.00 X 10®® 

i^^GC-H3 

4096® 

140.0 

3.44 X 10® 

1.37 X 10® 

44,679,543(« = 4) 

1.15 X 10®®(2 = 4) 


After the collapse of a DM halo, the hot gas gradually cools 
via radiative cooling, forming a cold gas disk at the halo center. 
Stars born from the condensed cold gas; and a stellar disk and a 
cold gas disk consist a galaxy (see section 2.4). The rate of gas 
cooling is calculated following the model proposed by White & 
Frenk (1991), which is adopted in most SA models. The time 
scale of radiative cooling, fcooi, is calculated as 






(3) 


2 prrip ni(r)A(rvir,^hot) ’ 
where ne(r) is the electron density of hot gas at r, Zhot is the 
metallicity of hot gas, and A is a metallicity-dependent cooling 
function provided by Sutherland & Dopita (1993). In each time 
step, the hot gas within the cooling radius cools and accretes 
onto the central cold gas disk. The cooling radius rcooi(l) is 
defined as the radius at which the cooling time scale is equal to 
the time elapsed since the halo formation epoch. If the cooling 
radius exceeds the virial radius Rvir, we set to rcooi = Rvu- In 


this case the mass accretion rate of cold gas should be limited by 
the free fall time, rather than the cooling time. However, we set 
the time step to be comparable to the dynamical time scale of 
halo at each epoch, thus this could not cause a serious problem. 

We further assumed that the radial profile of hot gas is kept 
unchanged until the DM halo doubles its mass, allowing the ex¬ 
istence of “cooling hole” at the centre of the halo (i.e., no hot gas 
is distributed at r < rcooi). This assumption is clearly unphys¬ 
ical, thus the effect on the gas cooling rate should be checked. 
Monaco et al. (2014) compared their SA model, MORGANA 
(Monaco et al. 2007), with other SA models and hydrodynam- 
ical simulations to test the cooling models. In their model, the 
cooling radius is treated as a dynamical variables and the gas 
profile is recomputed in each time step. They also consider the 
pressure balance between the hot gas and cooled gas, which de¬ 
termines the size of the cooling hole. One of other SA models 
examined in Monaco et al. (2014) adopts the cooling model of 
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White & Frenk (1991), as well as our model. Monaco et al. 
(2014) show that the different cooling models adopted in SA 
models only makes a marginal difference in cooling rate. See 
also De Lucia et al. (2010) for a test of the cooling models. 

As shown in equation (3), the cooling time scale depend on 
both the temperature and metallicity of the gas. In our model, 
the chemical enrichment of the hot gas due to the star formation 
and SN feedback is consistently solved as shown in subsection 
2.4. 

Note that the above assumption that the hot gas is heated 
up by shock at collapse of host halos is adopted just for sim¬ 
plicity. In reality, the cooling time scale of hot gas within 
galactic scale halos is much shorter than their dynamical time 
scale. Therefore, the hot gas should cool immediately rather 
than spherically re-distributing throughout the host halos. In 
any case, because the cooling time scale is very short, almost all 
the hot gas cools and thus our assumption is expected to work 
well. For its opposite case, within cluster-scale halos, the cool¬ 
ing time scale is very long owing to the high virial temperature 
and the AGN feedback. Again the assumption should be good. 
For the intermediate scale, we might need more sophisticated 
treatment. Along with the AGN feedback, these process should 
be improved in future versions of the model. 


2.3 Photoionization heating due to an UV 
background 

Intergalactic gas is photo-heated by the cosmological UV ra¬ 
diation field produced by galaxies and quasars. Because the 
heated gas cannot be accreted by small halos with shallow grav¬ 
itational potential wells, photo-heating quenches star formation 
in small galaxies and hence decreases the number densities of 
dwarf galaxies (e.g., Doroshkevich et al. 1967; Couchman & 
Rees 1986). The characteristic halo mass Me, below which a 
halo cannot retain the heated gas, has been investigated by using 
cosmological hydrodynamic simulations (e.g., Gnedin 2000). 

In this context, Okamoto et al. (2008) performed high- 
resolution cosmological hydrodynamical simulations with a 
time-dependent UV background radiation field. They found that 
the redshift evolution of the characteristic mass, Mc{z), is de¬ 
termined by the following factors for each halo: the relation 
between Tvir and the equilibrium temperature for the gas Teq at 
the edge of the halo, at which the density can be approximated 
as one third of the cosmic mean, and its merging history (see 
section 4 in Okamoto et al. 2008). They also found that the mass 
fraction of baryonic matter in halos with mass Mh at redshift z 
is well-fitted by the following formula, originally proposed by 
Gnedin (2000): 

/b(Mh,z) = (/b) 


X 


1 -h (2“uv 


/3 


1 ) 


Mb 
Me (2) 


— Q!UV 


3 

«UV 


(4) 


where the parameter auv controls the rate of decrease of fh 
in low-mass halos, here set to auv = 2. While /b(Mh, 2 ) 
equals to (/b) for the halos with Mb » Me (2), it goes to zero 
in proportion to (Mb/Me)® for the halos with Mb <C Me( 2 ). 
This decrease is attributed to the suppressed accretion of photo- 
heated baryonic matter onto the halos. This prescription, given 
by Okamoto et al. (2008), is newly incorporated into our u^GC 
model. Although all the above factors that determine Me ( 2 ) are 
evaluable in our u^GC model, we simply adopt their resultant 
Mc{z) itself in order to avoid a relatively large computational 
cost to obtain req((/b}pvir/3). 

The details how to incorporate the prescription of Okamoto 
et al. (2008) are as follows. Before reionization, which is as¬ 
sumed to instantaneously occur at 2 = 2reion, all halos contain 
baryonic matter with a mass fraction of fh = {fh) regardless of 
their masses, as described in subsection 2.2. After reionization, 
the expected baryon fraction fh of each halo with mass Mb that 
collapsed at 2 , denoted /b(Mh, 2 ), is calculated by equation (4) 
using a fitting formula for Me( 2 ): 


Mc(2) =6.5 X 10® 

X exp(—0.6042)exp [—( 2 / 8 . 37 )^^ ®] h~^ Mq. (5) 


This fitting formula is derived from the result of simulation 
of Okamoto et al. (2008) in which the reionization occurs at 
2 = 9.0. Figure 2 shows Me as a function of redshift. While the 
minimum halo mass of the lower resolution models (Mmin = 
8.79 X 10® h~^ Mq) is larger than Me, those of the higher res¬ 
olution models (Mmin = 1.10 x 10® h~^ Mq for the i/^GC-Hl 
and Mmin = 1.37 x 10® Mq for the ly’^GC-Hl, -H3, re¬ 
spectively) are smaller than Me at low redshift. In the halos 
with mass below Me, the gas heating by the UV background af¬ 
fects the properties of galaxies. To mimic this effect, Somerville 
(2002) also adopt the formulation of Gnedin (2000); however, 
they assume that Me (2) is given by constant Kire = 50 km s“^. 
In figure 2, we also show the redshift evolution of the halo mass 
with the fixed circular velocity of Ueirc = 17, 30 and 50 km s“^, 
respectively. We can see that the behavior of the Me proposed 
by Okamoto et al. (2008) and the Me given by fixed Kire are 
significantly different. 

In this paper we treat the effect of the gas heating due to 
the cosmic UV background as follows. For a halo with total 
baryonic mass (i.e., the sum of the masses of the stars, SMBH, 
cold gas, and hot gas of all galaxies in the halo) of Mb, tot < 
/b(Mh,2)M, the baryonic mass of fh(Mh,z)Mh — Mb, tot is 
added to the halo as hot gas with temperature of T^ir. On the 
other hand, for the halos with Mb,tot > fh{Mh,z)Mh, an ap¬ 
propriate mass of hot gas is removed from the halo keeping its 
metallicity unchanged so that the mass fraction of baryon in the 
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Fig. 2. The redshift evolution of the characteristic mass Me, beiow which a 
halo cannot retain the intergalactic gas due to the heating by the cosmolog¬ 
ical UV radiation field. The black filled circles show the results of the cos¬ 
mological hydrodynamical simulations performed by Okamoto et al. (2008). 
The thick solid line presents the fitting formula described in equation (5). The 
horizontal dashed lines denote the minimum halo mass of our fV-body sim¬ 
ulations (Mmin = 8.79 X 10®, 1.10 X 10® and 3.44 x 10® Mq from top to 
bottom, respectively). The thin solid lines correspond to the fixed circular 
velocity Vcirc = 50, 30 and 17 km s~^ from top to bottom. 


halo coincides with fh{Mh,z); this prescription mimics photo¬ 
evaporation by UV background radiation during the reioniza¬ 
tion. When Mb,tot — Mhot > /b(Mh, z)Mh, we have to reduce 
the additional cold gas masses in order to the mass fraction of 
baryon in the halo coincides with fh{Mh,z). However cold gas 
is much denser than hot gas, and may be self-shielded from the 
UV background radiation. Therefore we assumed that the cold 
gas component is not affected by the UV heating and allowed 
such halos to have larger baryon mass than /b(Mh, 2 )Mb. Note 
that the fraction of such halos is less than 1%, thus the treatment 
of the cold gas in this process does not have significant effects 
on the results presented in this paper. 

Although Okamoto et al. (2008) found that such photoevap¬ 
oration is important particularly just after the reionization (see 
middle panel of figure 5 of Okamoto et al. 2008), the effect in 
our model is assumed to occur for such less-massive halos with 
Mb.tot > /b(Mb, 2 )Mh, regardless of their collapsing redshifts 

when 2 < 2reion. 

2.4 Star formation and feedback in disk 

In this section we describe the star formation in cold gas disk 
and the reheating of cold gas by SNe. Our implementation fol¬ 
lows the standard recipe adopted in other SA models (e.g.. Cole 
et al. 2000). 

The cooling process of diffuse hot gas is followed by star 
formation in the cold gas disk. The star formation rate (SFR), tp, 
is given hy ip = Mcoid/Tstar, where Mcoia is the cold gas mass, 
and Tstar is the time scale of the star formation (SF). We assume 
that the star formation activity in galaxy disk is related to the 


dynamical time scale of the disk, ra = ra/Va, where ra and 
Va are the disk radius and disk rotation velocity, respectively, 
defined in subsection 2.8. Thus we adopt the following formula 
for star formation time scale Tstar, 


’7”star 


-1 

^star^^d 


i-f 


/ Vd 

VUhot/ 


( 6 ) 


where Sstar, ctstar, and Vhot are free parameters. Although the 
above modeling of the SF time scale well reproduces the sev¬ 
eral physical properties of observed galaxies as shown later (see 
section 5), there could be other modeling for the SF time scale. 
For example, we have examined the model in which the SF time 
scale depends on the global dust surface density of galaxy, and 
found that the choice of SF time scale model could have sig¬ 
nificant effect on the galaxy formation history (Makiya et al. 
2014). Although this would be promising for reproducing many 
aspects of observed galaxies, in this paper, we adopt this kind 
of standard prescription of star formation for simplicity. 

Consequent to supernova explosion, we assume that a frac¬ 
tion of the cold gas is reheated and ejected from galaxy at a rate 
of Mcoid/Treheat, whetc the time scale of reheating, Treheat is 
given as follows: 


"^stcir /^\ 

'T'reheat — ^ ; (, O 

with 

. , 8 , 

where Vkot and ahot are free parameters. 

With the above equations, we obtain the masses of hot gas, 
cold gas, and disk stars as functions of time (or redshift). The 
chemical enrichment associated with star formation and SN 
feedback is treated by extending the work of Maeder (1992). 
For simplicity, instantaneous recycling is assumed for SNe II, 
and any contribution from SNe la is neglected. 

In summary, the baryon evolution during the star formation 
process is described by the following equations: 


Mstar = aip{t), 

(9) 

Mhot = I3ip{t), 

(10) 

Mbh = fBnipit), 

(11) 

Mcoid = —{a + l3 + fBB)lpit), 

(12) 

(Mcold-^cold)' = [p — (a + /3 + /BH)2cold]'l/i 

(13) 

(Mhot^hot)'= PZcoldIp, 

(14) 


where Mstar and Mhot are the masses of stars and hot gas, 
respectively; ip = Mcoia/Tstar is SFR, Z^^id and Zhot are the 
metallicities of cold and hot gases, respectively, and Mbh is the 
mass of the nuclear SMBH. The constant parameter /bh con¬ 
trols the accretion rate of cold gas onto the SMBH during the 
starburst. In ordinary star formation process in disk, we assume 
that no cold gas gets accreted by the SMBH (i.e., /bh = 0.0). 
The galaxy merger and SMBH evolution will be detailed in later 














Publications of the Astronomicai Society of Japan, (2014), Vol. 00, No. 0 


7 


subsections (2.5 and 2.6). The locked-up mass fraction a and 
chemical yield p are chosen to he consistent with initial mass 
function (IMF). For the Chabrier IMF (Chabrier 2003), which 
is adopted in our standard model, a = 0.52 and p= 1.68 Z0. 

We can solve these equations analytically as 


AMcoid(f) 

AMstar(f) 

AMhot(^) 

^cold (f) 


M°m 


a + P 

P 

ct + p 
■^cold + P 


11 — exp 

AMcoid(^) 

AMcoid(^) 


-{a + P)- 


'^star 


^star 


MLtZLt + I 
a + P I 

P + Zcold{t)^ AMcold(f) 

— (■^cold(i) — ^cold).^cold Ij /AThot(^), 


(15) 

(16) 

(17) 

(18) 


(19) 


where the A symbol indicates that the variable is incremented 
or decremented in the current time step. All A variables are 
defined as positive. The superscript 0 denotes an initial values 
at the beginning of the time-step (i.e., t = 0). Note that here 
we assumed /bh = 0. For the case of burst-like star formation 
induced by major merger, see subsection 2.5. 


2.5 Mergers of galaxies and formation of spheroids 


After the merging of DM halos, the newly formed halo should 
contain two or more galaxies. The central galaxy in the most 
massive progenitor halo is designated as the central galaxy 
of newly formed halo, while the others are regarded as satel¬ 
lite galaxies. These satellite galaxies will fall into the central 
galaxy by dynamical friction (central-satellite merger). We set 
the merger time scale due to the dynamical friction as Tmrg = 
/mrgTfric, whcrc /mrg ~ 1 is adjustable parameter and rfric is the 
time scale of dynamical friction. If Tmrg is shorter than the time 
elapsed since the satellite galaxy enters the common halo, the 
satellite and central galaxy are merged. We reset this elapsed 
time to zero when the host halo mass doubles. 

For the time scale of dynamical friction, Tfric, we adopt the 
formulation of liang et al. (2008, 2010) which is obtained by 
fitting to the cosmological A'-body simulations : 


/(s) _ VciTcRcirc _ 

2C GMs In (1-fMh/Ms) 


( 20 ) 


where C = 0.43 is a constant fitting parameter, Kirc is the cir¬ 
cular velocity of the common halo, i?circ is the radius of the 
circular orbit of the satellite halo, and Mh and Ms are the total 
mass of the host halo and satellite halo, respectively. We sim¬ 
ply assumed that Rduc = Rh where i?h is the virial radius of 
the host halo. The function f{e) = O.OOe® "*^ -I- 0.60 accounts 
for the dependence of rfric on the orbital circularity e. We set 


to e = 0.5, which is the average value of e estimated by high- 
resolution A-body simulations performed by Wetzel (2011). In 
our A-body simulations we can resolve the satellite halos even 
after they entered the common halo, and therefore the time scale 
of dynamical friction can be directly drown from simulations; 
however, in this paper we adopted simplified formula described 
above fo save fhe compufational lime. The effect of this simpli¬ 
fication will be examined in a future paper. 

The satellite galaxies can randomly collide and merge 
(satellite-satellite merger). Makino & Hut (1997) conducted an 
A-body simulation of a system of same mass galaxies. They 
find that a merger rate, Icmh, is described by following simple 
scaling in this situation: 


fclVIH = 


^ 1 

/lMpc\3 1 

^ ^gal \ 

500 ' 

[ Rh J ' 

1 0.1 Mpc ) 


/ '^gal \ ** ( 

VlOOkms-1/ 


300 km s 

^halo 


3 

Gyr-\ 


( 21 ) 


where A, o-gai, Tgai and crhaio are the number of satellite galax¬ 
ies, one-dimensional velocity dispersions of the galaxy, galaxy 
radius and parent halo, respectively. In our model a satellite 
galaxy will collide with another satellite galaxy picked out at 
random, with the probability of At x Icmh, where Af is the time 
step of the calculation. 

We consider two distinct modes for galaxy merger, i.e., ma¬ 
jor merger and minor merger. If the ratio of baryonic mass 
(stars, cold gas, hot gas and SMBH mass) of two merging galax¬ 
ies, / (< 1), exceeds the critical value /bulge, major merger oc¬ 
curs. Major mergers induce burst-like star formations, in which 
all of the cold gas in the merging system turns into stars and hot 
gas. The star formation and SN feedback law is the same with 
the disk star formation (see section 2.4), except for assuming the 
very short star formation time scale (rstar —>■ 0). The bulges and 
stellar disks of the progenitor completely reform into the bulge 
component of the new galaxy, together with the stars born dur¬ 
ing the merger. Note that when applying the SN feedback law, 
the disk velocity Vd is replaced by the velocity dispersion of the 
new bulge 14 (defined in subsection 2.8). 

On the other hand, if / < /bulge, a minor merger occurs. 
In this case, stellar and cold gas components of the smaller 
galaxy are absorbed into the bulge and cold gas disk of the larger 
galaxy, respectively, with no starburst events. 


2.6 Supermassive black holes 

Along with the evolution of galaxies, SMBHs at galaxy cen¬ 
ters also evolve by the following mechanismas: (1) SMBH co¬ 
alescence, (2) accretion of cold gas (during a major merger of 
galaxies), and (3) “radio-mode” gas accretion. Note that we as¬ 
sume a central SMBH in every galaxy. When the galaxies first 
collapsed, the seed BH have formed with mass Mgeed, which is 
a tunable parameter. 
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It has been shown by theoretical studies that a major 
merger of galaxies can drive substantial gaseous inflows into 
a galaxy center (e.g., Mihos & Hernquist 1994, 1996; Barnes 
& Hernquist 1996; Di Matteo et al. 2005; Hopkins et al. 2005, 
2006). We assume that a fraction of this inflowing cold gas gets 
accreted by the central SMBH. The mass of cold gas accreted 
by the SMBH, AMbh is modeled as follows: 

AMbH = /BnAMstar,burst, (22) 

where /bh is a constant, and AMstar,burst is the total mass 
of stars formed during the starburst. We set /bh = 0.005 to 
match the observed relation between masses of host bulges and 
SMBHs at 2 : = 0 (see subsection 5.4). The accretion of cold 
gas triggers the quasar activities. For more detailed model 
of quasars, see Enoki et al. (2003), Enoki et al. (2014), and 
Shirakata et al. (2015). 


2.7 AGN feedback 

To reproduce the observed break at the bright end of the lumi¬ 
nosity functions (LFs), we introduced the so-called radio-mode 
AGN feedback process into our model. In this mode, the hot gas 
accreted by SMBH powers radio jet that injects energy into the 
hot halo gas, quenching the cooling of the hot gas and resultant 
star formation in the massive halo. Radio-mode AGN feedback 
is also expected to contribute to the downsizing evolution of 
galaxies (e.g., Croton et al. 2006; Bower et al. 2006). 

Our implementation of AGN feedback follows the formula¬ 
tion of Bower et al. (2006). In their formulation, gas cooling in 
the halo is inhibited when the following conditions are satisfied: 

Q^coolldynfrcool) ^ fcool (27) 

and 

eSMBHT/Edd > L cool ? (28) 


Considering the very short time scale of starburst (f/rstar —^ 
oo) assumed here and the mass accretion onto the nuclear 
SMBH, we solve equations (9)-(14) to obtain the following: 


AMatar = 

AMhot = 
AMbh = 

A(Mhot^hot) = 


a + j3 + /bh 

P 

a + P + /bh 


M'cold, 


M?old, 


/bh 


a + P + /bh 

P 

a + P + /bh 
P 

01 + P -{■ /bh 


M'cold, 


+ -^cold ) Mcold, 


(23) 

(24) 

(25) 

(26) 


where AMstar, AMhot, AMbh, and A(Mhot^hot) are the in¬ 
creasing amount of the stellar mass, hot gas mass, BH mass and 
the metal mass in the hot gas, respectively, during a starburst. 
The superscript 0 indicates the total values in the merger pro¬ 
genitors. We again emphasize that all the cold gas is exhausted 
in our starburst model. 

During a merger event, an SMBH also increases its mass via 
the SMBH-SMBH coalescence. In this paper, we simply as¬ 
sume that SMBHs merge instantaneously right after the merger 
of their host galaxies, because it is difficult to estimate the time 
scale of SMBH mergers owing to the existence of many compli¬ 
cated physical processes such as the dynamical friction, stellar 
distribution, multiple SMBH interaction, and gas dynamical ef¬ 
fects (see, e.g., Colpi 2014). As shown in Enoki et al. (2004), 
the mass growth of SMBHs in our model is mainly due to the 
gas accretion during major merger, at least, at z < 1, and there¬ 
fore therefore the assumption of the instantaneous coalescence 
does not have significant effects. The other evolution channel, 
radio-mode gas accretion related to the AGN feedback process, 
is described next. 


where fdyn is the dynamical time scale of the halo at the cool¬ 
ing radius, fcooi is the time scale of gas cooling, Z/Bdd is the 
Eddington luminosity of the AGN, Lcooi is the cooling luminos¬ 
ity of the gas, and Ofcooi and ssmbh are free parameters that are 
tuned to reproduce the observations. Under these conditions, 
AGN feedback is limited to haloes in quasi-hydrostatic equi¬ 
librium, and having a sufficiently evolved SMBH. In the halo 
experiencing AGN feedback, the SMBH at the center grows by 
accreting hot halo gas. Bower et al. (2006) assumed that the 
accretion flow is automatically adjusted by itself so that heat¬ 
ing luminosity balances with cooling luminosity, namely, the 
accretion rate is set to Mbh = Tcooi/??c^. Here 77 is the radia¬ 
tive efficiency. We assumed 77 = 0.1 for all SMBHs, based on 
the observational estimation of Davis & Laor (2011). The value 
of 77 does not significantly affect on the results since the mass 
growth of SMBHs is dominated by the gas accretion during ma¬ 
jor merger. 

2.8 Size of galaxies and dynamical response to gas 
removal 

This subsection explains how we estimate the galaxy size, disk 
rotation velocity and velocity dispersion of bulges. Our recipe 
of size estimation almost follows the procedure of Cole et al. 
( 2000 ). 

2.8.1 Disk formation from cooled gas 

First, we estimate the size of galaxy disk as follows. We assume 
that the hot halo gas has the same specific angular momentum as 
the DM halo and collapses to the cold gas disk while conserving 
the angular momentum. We introduce the dimensionless spin 
parameter Ah as Ah = L\E\^^‘^, where L is the angu¬ 
lar momentum, E is the binding energy and M is the DM halo 
mass. Although the distribution of Ah is often approximated 
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by a log-normal distribution (e.g., Mo et al. 1998), it has been 
known that the distrihution of Ah deviates from log-normal in 
large M-hody simulations (e.g., Bett et al. 2007). However, ac¬ 
cording to Bett et al. (2007), the shape of distrihution depends 
on halo finding algorithm and the log-normal function is still 
slightly better for FoF halos than the modified fitting function 
proposed by them. In this paper, we simply adopt the log¬ 
normal distrihution 


p(AH)dAH = —j =—exp 
vSttcta 


(In Ah - In A)^ 
2ol 


dlnAn, 


(29) 


where A and a\ denote the mean and logarithmic variance of 
the spin parameter, respectively. Here we use A = 0.042 and 
a\ = 0.26., which are the value obtained by Bett et al. 2007 for 
FoF halos. 

Using the spin parameter Ah, the effective radius of a 
resultant cold gas disk is expressed as follows (Fall 1979; Fall 
& Efstathiou 1980; Fall 1983): 


Ti = (1.68/V^)AHi?i, 


(30) 


where the initial radius of the hot gas sphere, Ri, is set to the 
virial radius of the host halo or the cooling radius, whichever 
is smaller. In each time step, the disk size of central galaxies 
are updated if their disk mass have increased from the previous 
time step. At this time, we set the disk rotation velocity Vd to 
be the circular velocity of its host halo. 


2.8.2 Dynamical response to disk star formation 

After the formation of rotationally supported disks, the SN feed¬ 
back subsequent to disk star formation expels cold gas continu¬ 
ously. As the baryonic mass of galaxies decreases, the grav¬ 
itational potential well becomes shallower, depending on the 
mass ratio of baryons to DM within the galactic disk. In re¬ 
sponse to the variation of the depth of the gravitational poten¬ 
tial well, gravitationally bound systems expand and their rota¬ 
tion speed slows down (Yoshii & Arimoto 1987). We refer to 
this effect as the dynamical response here. Dwarf galaxies hav¬ 
ing shallow gravitational potential wells and therefore suffered 
significant SN feedback are affected more by the dynamical re¬ 
sponse. Using our SA models taking into account this for star- 
burst, we have shown that this affects the scaling relations of 
elliptical galaxies especially for dwarfs (Nagashima & Yoshii 
2004; Nagashima et al. 2005). See those papers for the scheme 
of introducing the dynamical response in SA models. In the 
present paper, we also apply this effect to disk evolution. 

The basic result for disks used here is given by Koyama et 
al. (2008). At first, we assume a galactic disk within a static 
DM halo and approximate the density distributions of disks and 
DM as the Kuzmin disk (Kuzmin 1952, 1956) and the Navarro- 
Frenk-White (NFW) profile (Navarro, Frenk & White 1997), 
respectively. Then, we consider that the gas mass of disks grad¬ 
ually decreases due to the SN feedback, that is, the so-called 


adiabatic mass-loss, and that the dynamical response to the gas 
removal on the disks. The initial radius of cold gas disk is cal¬ 
culated by equation (30). Here we assume that the disk size is 
determined only by the gravitational potential of the host dark 
halo, conserving the specific angular momentum of the cooled 
gas. Although this might be too simple because the central re¬ 
gion of galaxies would form dynamically with cooling gas, it 
should be a good approximation for the outer disk. Thus we 
take this treatment in this paper as usual. 

Here we define A4 and R as ratios of mass and size at a final 
state relative to those at an initial state, and Zi and z/ as ratios of 
baryonic disk size relative to size of dark halos at those states. 
According to Koyama et al. (2008), we obtain 

= l + (31) 

R rrii 

where mi is the mass ratio of baryons to dark matter at the ini¬ 
tial state, and q{z) is a function depending on the distributions 
of baryons and dark matter. In this case, we cannot obtain an 
analytic form of q{z). Instead, we expand the above equation 
around 2 = 0 and i? — 1 ~ 0 as 


A1 = -+D(7?-l), 
where 


(32) 


D =— 

mi L 


In(l-l-c) — 


1 + c. 

16 2 3 34 

-C Zi — C Zi 

3 


^3-|-21n^ 

V 8 2 2 / 


(33) 


and c is the concentration parameter described in section 2.2. 
Note that we take a higher order term of Zi for q{z) than that 
written in equation (A5) in Koyama et al. (2008). 

The approximation used here is justified as follows. It is 
expected that the change in sizes and disk rotation velocities 
during a time-step is very small because of the quiescent star 
formation, and that the size of baryonic disks is smaller suf¬ 
ficiently than that of dark halos. These mean i? — 1 <C 1 and 
Zi "C 1. We have checked that these assumptions are indeed 
validated in our model. 

The change of the disk rotation velocity is given by 

1/2 

(34) 




mf/zf+4:f{zf) 


mi/zi + 4:f{zi) 


where f{z) is also a function depending on the distributions of 
baryons and dark matter, similar to q{z). The form of f(z) is 
shown in equation (A4) in Koyama et al. (2008). 

We would like to recall here the well-known results for non¬ 
dark matter case, R = V~^ = A4~^. These relations are ob¬ 
tained by setting f{z) and q{z) to zero. In the opposite limit¬ 
ing case, because dark matter dominates, q{z) becomes much 
larger. In this case, even if A4 becomes zero, R and U do not 
vary. This corresponds to the case discussed in Dekel & Silk 
(1986). 

The effect of the dynamical response on the disk shall be 
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discussed in detail in another paper. 


2.8.3 Dynamical Response to starburst and spheroidal rem¬ 
nants 

The size of the bulge formed in a major merger is characterized 
by the virial radius of the baryonic component. Applying the 
virial theorem the total energy in each galaxy is calculated as 
follows: 


Ei = + {Md,i + Mcoid)K?,i], (35) 


where Mb, Mbh, Md and Mcoid are the masses of the bulge, 
central black hole, stellar disk and cold gas disk, respectively, 
and Vb and Vd denote the velocity dispersion of the bulge and 
the rotation velocity of the disk, respectively. The subscript 
i = {0,1,2} indicate the merged galaxy, larger progenitor galaxy 
and smaller progenitor galaxy, respectively. Furthermore the or¬ 
bital energy Emh between the progenitors just before the merger 
is given as follows: 


F/orb — 


E 1 E 2 

{M2/Mi)Ei + {Mi/M2)E2 ■ 


(36) 


By energy conservation, we obtain the following: 


fdiss {El + E 2 + Eorh) — Eo, 


(37) 


where fdiss is the fraction of energy dissipated from the system 
during major merger. The rate of energy dissipation depends on 
complicated physical processes such as the viscosity and fric¬ 
tion due to gas. In this paper, we simply parameterize fdiss as 
follows: 


fdiss — 1 “t“ t^diss/g; 


(38) 


gas removed from galaxies and the mass of the system changes. 
At this time the structural parameters of galaxies also changes 
due to the dynamical response. We include this effect into our 
model adopting the Jaffe model (Jaffe 1983). In this paper 
we assume the case of slow (adiabatic) gas removal compared 
with the dynamical time scale of the system, similar to that for 
disks. For the case of rapid gas removal, we refer the reader 
to Nagashima & Yoshii (2003). According to the Nagashima 
& Yoshii (2003), the effect of dynamical response becomes 
stronger in the case of non-adiabatic gas removal. Therefore 
the assumption of the adiabatic gas removal should be consid¬ 
ered as conservative. We should keep in mind that the effect 
might be stronger for dwarf ellipticals having shorter time scale 
of gas removal compared to giants. 

Defining hy M,R and U the ratios of mass, size and velocity 
dispersion at a final state relative to those at an initial state, the 
response under the above assumption is approximately given by 


Tf l + D/2 
~ n MED 12' 

U=Yhl^ l MIR + Df{zf )72 
Vb,i V l + Df{zi)/2 


(41) 

(42) 


where D = l/yiZ^,y and z are the ratios of density and size of 
baryonic matter to those of dark matter. We use equation (36) 
in Nagashima et al. (2005) for the form of f{z). The subscripts 
i and / stand for the initial and final slates in the mass loss pro¬ 
cess. Note that U is the ratio of velocity dispersion, different 
from that for disks. The contribution of dark matter is estimated 
from the central circular velocity of halos, Kent, which is de¬ 
fined below. 


where 


r _ Mcoid 

Mstar + Moold + MbB 

is Ihe gas mass fraction of the merging system and Kaiss is a di¬ 
mensionless parameter. Flere we set Kaiss = 2.0 to reproduce the 
distribution of size and velocity dispersion of elliptical galaxies 
(see section 5.2). There are several studies on this issue by using 
hydrodynamical simulations and SA models, and it is confirmed 
that above parameterization of fdiss can be a good approxima¬ 
tion (see, e.g., Hopkins et al. 2009; Shankar et al. 2013). 

We assumed that there remains only the bulge compo¬ 
nent supported by velocity dispersion just after the merger. 
Therefore the velocity dispersion and the size of merger rem¬ 
nant can be estimated from following equations, 

F;o = -iMtot.oK'o, (39) 

and 


GMtot,o 


(40) 


where Mtot,o is the total baryonic mass of the merger remnant. 
As a consequence of star formation and SN feedback, part of 


2.8.4 Back reaction of dynamical response to dark halos 

When galaxies suffer the dynamical response to gas removal 
caused by the SN feedback, the dark matter within a central 
region of dark halos hosting the galaxies must also suffer the 
dynamical response as its back reaction. For simplicity we com¬ 
pute the dynamical response on the dark matter distribution af¬ 
ter the computation of the dynamical response on baryons, al¬ 
though they occur simultaneously in reality. Here we ignore 
the effect of the so-called adiabatic contraction for dark mat¬ 
ter during the condensation of cooled gas. This is because the 
central region of galaxies should form not adiabatically but dy¬ 
namically. Thus we assume that the cooled gas condenses and 
relaxes dynamically together with dark matter and is removed 
adiabatically by the SN feedback affecting the central region of 
dark halos as the back reaction. This would require detailed 
research by using hydrodynamic numerical simulations. 

Here we focus on the region within the half-mass radius of 
central galaxies, at which the density of baryons is expected to 
be comparable to that of dark matter. To take into account this 
process, we define a central circular velocity of dark halo Kent, 
approximately within the effective radius of the central galaxy. 
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When a dark halo collapses without any progenitors, Kent is set 
to Kirc After that, although the mass of the dark halo grows 
by subsequent accretion and/or mergers, Kent remains constant 
or decreases with the dynamical response. When the mass is 
doubled, Kent is set to Kirc again. According to Nagashima 
& Yoshii (2004) and Nagashima et al. (2005), we assume that 
Kent is lowered by the dynamical response to mass loss from a 
central galaxy of a dark halo by SN feedback as follows, 

Kent,/ _ Mf /2 + Md{ri/rd) , 

Kent,i “ Mi/2 + Md(ri/rd)' ^ ’ 

The change of Kent in each time-step is only a few per cent. 
Under these conditions, the approximation of static gravita¬ 
tional potential of dark matter is valid even during starbursts. 
This also applies to subhalos. Rigorously speaking, we must 
assume an isothermal distribution of dark matter, in which the 
density is proportional to the inverse of r^, because the above 
equation indicates the dynamical response to mass loss within 
the half-mass radius of central galaxies. In spite of this, this 
should be good approximation because the NFW profile has a 
slope —1 and —3 within and outside the core radius, respec¬ 
tively, which means that we can expect that the effective slope 
would be approximately —2. Of course this expectation is op¬ 
timistic since we considered the inner region of a halo where a 
slope is —1. However, we need detailed hydrodynamical simu¬ 
lations to know the actual mass profile since the adiabatic con¬ 
traction due to gas cooling would affect the slope. 

Once a dark halo falls into its host dark halo, it is regarded 
as a subhalo. Because subhalos do not grow in mass in our 
model, the central circular velocity of the subhalos monotoni- 
cally decreases. Although the change of Kent during a time- 
step is small, accumulated change cannot be negligible owing 
to the monotonicity. Therefore, this affects the time scales of 
mergers. 

The details of the dynamical response are shown in 
Nagashima & Yoshii (2003, 2004) and Nagashima et al. (2005) 
for bulges and Koyama et al. (2008) for disks. The effect of 
the dynamical response is the most prominent for dwarf galax¬ 
ies of low circular velocity because of the substantial removal 
of gas due to strong SN feedback (Yoshii & Arimoto 1987; 
Nagashima & Yoshii 2004). If the dynamical response had not 
been taken into account, velocity dispersions of dwarf ellipticals 
would have been much larger than those of observations, deter¬ 
mined only by circular velocities of small dark halos in which 
dwarf ellipticals resided. For giant ellipticals, on the other hand, 
the effect of the dynamical response is negligible because only a 
small fraction of gas can be expelled due to weak SN feedback. 
Similarly, for disks, in order to reproduce the observed Tully- 
Fisher relation, the dynamical response on disks is required. 
Otherwise, the slope for dwarf spirals becomes different from 
observed one as shown in Nagashima et al. (2005). This point 
will be discussed in detail in another paper. 


2.9 Photometric properties and morphoiogicai 
identification 

Calculating the baryonic processes described in the above sub¬ 
sections, we finally obtain the SF hand metal enrichment histo¬ 
ries of each galaxies. From this information, we can calculate 
spectral energy distribution (SED) of model galaxies by using a 
stellar population synthesis code of Bruzual & Chariot (2003). 

To estimate the extinction of starlight, we first assume that 
the dust-to-cold gas ratio is proportional to the metallicity of 
the cold gas; second, we assume that the dust optical depth is 
proportional to the dust column density. The dust optical depth 
Tdust is then calculated as follows: ruust is given by 

—<«' 

where ra is the effective radius of the galaxy disk and tq is a tun¬ 
able parameter that should be chosen to fit the local observations 
(such as LFs). The wavelength dependence of optical depth is 
assumed to follow Calzetti-law (Calzetti et al. 2000). Dust dis¬ 
tribution is assumed to obey the slab dust model (Disney et al. 
1989) for disks. 

In our model, a major merger induces starburst activity, in 
which all the cold gas turns into stars and hot gas. Therefore, no 
cold gas and dust exist immediately after the starburst. Hence, 
the dust optical depth exactly equals to zero and galaxy color 
becomes too blue compared to the observations. To avoid this 
problem, we estimate the amount of dust extinction during the 
starburst as follows. First, we randomly assign the merger 
epoch within the current time step. Second, we calculate the 
amount of remaining dust at the end of the time step. At this 
time, the time scale of gas consumption during the burst is 
assumed as the dynamical time scale of the merged system, 
rh/Vb- The dust geometry is assumed to be the screen model. 

The morphological types of model galaxies were determined 
by the bulge-to-total luminosity ratio in the B-band. In this pa¬ 
per we follow the criteria of Simien & de Vaucouleurs (1986): 
galaxies with B/T > 0.6, 0.4 < B/T < 0.6, and B/T < 0.4 are 
classified as elliptical, lenticular, and spiral galaxies, respec¬ 
tively. According to Kauffmann & White (1993) and Baugh et 
al. (1996), this classification reproduces well the observed type 
mix. 

3 Parameter settings 

As described in section 2, our model is constructed from physi¬ 
cally motivated prescriptions of several astrophysical processes. 
However, a number of free parameter remain. Here we describe 
the parameter setting procedure. 
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3.1 Overview of parameter settings 


3.2 Markov chain Monte Cario anaiysis 


For the cosmological parameters, we adopt the Planck cosmol¬ 
ogy (Planck Collaboration et al. 2014). The several free pa¬ 
rameters related to astrophysical processes are listed in table 
2. Seven of these parameters, namely, Ostar, Tstar, Ohot, Vhot 
Qcooi, esMBH and Mseed were tuned to fit the local optical (r- 
band) and near IR (if-band) LFs and the local mass function 
(MF) of cold neutral hydrogen, by using a MCMC method (see 
next section). We use the local LFs and Fit MF as the fiducial 
references in model calibration, since they are robustly deter¬ 
mined from recent large and deep surveys. The other parame¬ 
ters, /bulge, /mrg, /bh , Tvo and Kdiss, atc manually tuned by 
comparing other observations, since they cannot be constrained 
by the local LFs and H i MF. 

The galaxy merger-related parameters, /bulge and /mrg, are 
closely related to the abundance of elliptical galaxies, hence 
they can be constrained by the LFs divided by morphological 
class. However there are still some uncertainties in the deter¬ 
mination of morphology, thus we did not use them in the fit¬ 
ting. In this paper we simply assumed that /bulge = 0.1 and 
/mrg = 0.8, which is the same value with NOS. The mass frac¬ 
tion accreted by SMBH during a starburst, /bh, affects the 
bright-end shape of LFs through the AGN feedback; however 
it is degenerated with other AGN feedback-related parameters, 
esMBH and Mseed, and is poorly constrained by LFs. Thus we 
tuned /bh to reproduce the observed BH mass - bulge mass 
relation and mass function of SMBHs, which are significantly 
affected by /bh but not by other two parameters. We have found 
that /bh = 0.005 is suitable to reproduce the observations in the 
case of /bulge = 0.1 and /mrg = 0.8 (see section 5.4). The co¬ 
efficient of dust extinction, tvo, was set to the value adopted in 
NOS, namely tvo = 2.5 x 10“®. The parameter related to the 
energy loss fraction in a major merger (Kdiss = 2.0) was chosen 
to fit the size-magnitude relation of elliptical galaxies (see sec¬ 
tion 5). Throughout this paper, we adopt the Chabrier IMF in 
the mass range 0.1-100 Mq. 

Although the MCMC method is numerically economical, it 
still requires approximately ~ 10® realizations to estimate the 
reliable parameter range. Therefore, to restrain the runtime of 
each realization within a few seconds, we employed the i/^GC- 
SS model for the A-body data in the MCMC fitting, which has 
the lowest mass resolution and the smallest box size. The mass 
resolution of A-body data could have complicated effects on 
the merging history of DM halos, thus there is no guarantee 
that the parameters tuned for u^GC-SS model work well for 
other A-body runs. However no significant differences between 
the i/^GC-SS and i^^GC-H2 (the highest resolution model) are 
found in the r- and A-band LFs and H i MF in the magnitude 
and mass range used in the fitting (see section 3.4 and figure 4, 

5). 


MCMC analysis was implemented by the Metropolis-Hastings 
algorithm (Metropolis et al. 1953; Hastings 1970), which is the 
most commonly used MCMC method. This method requires 
the proposal distribution g, which suggests a candidate point 
for the next step, given the previous sampling point. We as¬ 
sume a Gaussian distribution function for q. The variance of 
the Gaussian function is manually selected, to decrease the con¬ 
vergence time. We run 8 MCMC chains in parallel from ran¬ 
dom starting points. Each chain has about 50,000 realizations, 
excluding initial 10,000 steps of ”burn-in” phase. The conver¬ 
gence of MCMC chain is checked by the Gelman-Rubin diag¬ 
nostic test (Gelman & Rubin 1992). In this method, the dif¬ 
ference between the multiple MCMC chains are quantified by 
the ratio of the variance between chains to the variance within 
chain, R. In this paper the chain is considered to have converged 
when R < 1.1. As a result of the fitting, all the free parameters 
examined here reach convergence. We simply assume the uni¬ 
form distribution for the prior probability distribution, with the 
range listed in table 2. Although the bounds of prior distribu¬ 
tions are physically chosen, the ranges are set to be wide in 
order to cover a large model space, since our knowledge about 
the posetrior distribution of parameters is limited. 


3.3 Observational data and error estimation 


In this subsection, we describe the observational data used in 
the MCMC fitting. The local r- and A-band LFs were obtained 
by the Galaxy and Mass Assembly (GAMA) survey (Driver et 
al. 2012), and H i MF was extracted from the data of the Arecibo 
Legacy Fast ALFA (ALFALFA) survey (Martin et al. 2010). 
For each realization, the likelihood is calculated as follows: 

E = £oexp(^-2^^±24±24^, (45) 

where Co is an arbitrary constant and Xr, Xk Xhi are the x^ 
values of the r-band LF, A-band LF and Hi MF, respectively. 
These values are estimated as follows: 


X^('^obs|6') 


\ '' {^i,obs 0i,model(^)} 

7. ’ 


(46) 


where (j)i,oha denotes the value in the ith bin of the observed 
LF (or Hi MF), </i,model(0) is the value of the model in the ith 
bin obtained with the parameter set 6 and (7i,obs and (Ti,model 
are the errors in the observation and model in each bin, respec¬ 
tively. The errors in the observed LFs only include Poisson er¬ 
rors (Driver et al. 2012), while the errors in the observed H i MF 
include the systematic errors in mass estimation in addition to 
Poisson errors. 

The errors in the model predictions, cri^modei, were assumed 
as the sum of Poisson statistical errors and systematic errors 
coming from cosmic variance. Although the most of SA mod- 
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els assume Poisson errors for the model (e.g., Henriques et al. 
2009; Lu et al. 2014), it is controversial whether this assump¬ 
tion is appropriate or not. However typical value of Poisson er¬ 
rors are less than 1 % of the error coming from the cosmic vari¬ 
ance described helow, it does not have significant effect on the 
parameter fitting. Although the errors of the observed H i MF 
include the systematic errors come mainly from uncertainty of 
mass estimation, especially for low H i mass galaxies (see, e.g., 
Zwaan et al. 2005; Martin et al. 2010), we does not include it in 
the errors of our model. 

The effect of cosmic variance is estimated as follows. First, 
we ran the model using the i/^GC-S for the A-body data, which 
has a larger box (L — 280h~^ Mpc) than that in the MCMC fit¬ 
ting (L = 70h~^ Mpc), and randomly picked out the L = 70h~^ 
Mpc box from the large box data in ~100 trials. Following this, 
we drew the LFs and H i MF from the small boxes and deter¬ 
mined their uncertainties (approximately 20%) in each bin. We 
accounted for this 20% uncertainty in (7i,moodei, in addition to 
Poisson errors. 

Populations of dwarf galaxies with low surface brightness 
are known to exist, and the faint-end slope of observed LFs 
may be affected from the surface brightness limits of galaxy 
surveys (e.g. Blanton et al. 2005). According to Baldry 
et al. (2012), the incompleteness of GAMA samples become 
larger than 30% at pr, 5 o ^ 23.5 mag arcsec”^, where pr,so 
is the surface brightness within the Petrosian half-light radius. 
Therefore, we adopt this limit in calculating the model LFs. To 
calculate the Petrosian surface brightness, we require the light 
profile of galaxies. However, because our model does not re¬ 
solve the internal structure of galaxies, we converted the effec¬ 
tive radius and total magnitude into the Petrosian radius and 
Petrosian magnitude, respectively, fixing the Sersic index ris of 
bulge (ris = 4) and disk (ns = 1) components for all galaxies. 

The mass of cold atomic hydrogen of model galaxies are 
estimated as follows. First we assume that the 75% of the cold 
gas is composed of hydrogen. This cold hydrogen will be split 
into atomic and molecular; however, our current model does 
not follow the complex history of the formation of molecular 
hydrogen. Therefore we simply assume a fixed H 2 -to-Hi ratio 
for all galaxies. According to the observational estimation of 
Keres et al. (2003) and Zwaan et al. (2005), a global mass ratio 
of molecular to atomic hydrogen is ~ 0.4. Thus the mass of 
cold atomic hydrogen is estimated as 

Mhi = 0.75/(l + 0.4)Mcoid. (47) 

Note that similar approach was used in other SA models (e.g. 
Power et al. 2010; Lu et al. 2012). When fitting Hi MF, we 
only used the data points acquired at Mhi > 10 ® Mq, because 
at masses below this limit, the mass resolution the A-body data 
would affect the shape of low mass end of H i MF (see next sub¬ 
section). The uncertainties in the observed Hi MF also increase 


below this limit due to the incompleteness of the survey (Martin 
et al. 2010 ). 

3.4 Fitting resuits 

The diagonal panels of figure 3 present the ID posterior proba¬ 
bility distributions of the parameters tuned in the MCMC fitting. 
From the ID posterior probability distributions, we computed 
the medians and 10 and 90 percentiles of each parameter; the 
statistics are summarized in table 2. The off-diagonal panels 
of figure 3 present the 2D posterior probability distributions of 
all combinations of the seven free model parameters (grey con¬ 
tours). The ID distributions of the five parameters, oistar, Tatar, 
ahot, Vhot, and Qcooi are highly-peaked, indicating that they are 
well constrained within the assumed range. On the other hand, 
esMBH and Mseed have broad distribution. This can be under¬ 
stood as follows. If these parameters are large enough, the sec¬ 
ond condition of AGN feedback [equation (28)] will be satisfied 
in all halos. In such case, the specific value of these parameters 
no longer affect the shape of LFs and MF, and therefore only 
the lower boundary of these can be constrained from the fitting 
of LFs and H i MF. The posterior distribution of Mseed suggests 
that the Mseed should be larger than 10 ®Mq. This seed BH 
mass is somewhat higher than other SA models. In our model 
a fraction of central galaxies is bulge-less (i.e. they have never 
experience any major merger event till z = 0). The SMBH mass 
is equal to Mseed in such galaxies, thus large value of Mseed is 
required in order to make AGN feedback work in such halos. 

Figure 4 present the r- and A-band LFs in the model with 
the MCMC-obtained best fit parameters. The model closely 
matches the observations over all magnitude ranges. The 
shaded regions indicate the la error in the model LFs, esti¬ 
mated from the la confidence interval of each parameter. To 
see the effect of the mass resolution, we also plot the results of 
i/^GC-H2 model with the same parameters. These models are 
consistent within the Icr error. 

Figure 5 shows H i MF computed by the best-fit model. Data 
below the lower limit of the H i mass (solid vertical line) were 
excluded in the MCMC fitting because they deviated when the 
model was run at higher mass resolution (u^GC-Hl model, 
Mmin = 1.37 X lO®h~^M0; dashed line). Although there re¬ 
mains uncertainties in both the model and the observation, the 
model seems to underpredict the abundance of lower H i mass 
galaxies (Mhi < 10 ® M©). The similar trend is seen in other 
SA models. For example, Gonzalez-Perez et al. (2014) find that 
their model also underpredicts the abundance of galaxies at the 
lower-mass end of Hi MF (see also Lagos et al. 2014). They 
conclude that this is mainly due to the limited mass resolution 
of their A-body data. However, even in the i^^GC-H2 model, 
which has approximately two orders of magnitude higher mass 
resolution than that of Gonzalez-Perez et al. (2014), the lower- 
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Fig. 5. H I mass function of the best fit modei. The biack solid line represents 
the model (i/^GC-SS) with the best-fit parameters determined by MCMC 
fitting. The black dotted line shows the model with the same parameters but 
using the high-resolution JV-body data (!^^GC-H2). The mass resolution of 
w-body data would affects the shape of MF below Mm ~ 10®Mq (shown by 
vertical solid line); therefore, data below this were excluded in the parameter 
fitting. The shaded region denotes the lo- error in the model, estimated from 
the posterior probability distribution of the parameters (figure 3). Black filled 
circles are the observational data obtained by ALFALFA survey (Martin et al. 
2010 ). 

mass end of the H i MF is still underpredicted. This result might 
suggest that the more realistic modeling of star formation and 
SN feedback is required (e.g., Lu et al. 2014; Benson 2014). 
Furthermore, non-virialized gas which is not included in our 
model and/or FI i gas with low H i column densities below the 
observation limits of the current Fit blind surveys might con¬ 
tribute to the low end of the HI MF (e.g., Okoshi et al. 2010). 
We will further investigate this issue in the future. 

4 Numerical galaxy catalog 

Following the above procedures, we finally obtained the numer¬ 
ical galaxy catalog. This catalog contains various data on each 
mock galaxies: redshift, three-dimensional positions, physical 
quantities such as stellar mass, gas mass, metallicity, star forma¬ 
tion rate, effective radius, and magnitudes in several passbands 
in the UV-NIR regime. More information is provided at here*. 

Figure 6 plots the spatial distribution of the model galax¬ 
ies from 2 = 0.0 to 2 = 11.6 (corresponding to approximately 
lO'* Mpc along the comoving radial distance), plotted on the 
past light-cone of an observer at 2 = 0. Here we show the result 
of the z 7^GC-H1 model. The light-cone is generated by patch¬ 
working the model outputs at various redshift slices. During 
the patchworking, the simulation box was randomly shifted and 
rotated to avoid artifacts in the spatial structure. Web-like struc¬ 
tures are clearly visible in this figure. Thanks to the high mass 
resolution of the model, we can observe galaxies even at 2 > 10. 


Fig. 7. Effective radius of spiral galaxies plotted against /-band magnitude. 
The blackfilled squares with error bars show the median and the 10th to the 
90th percentile of the predicted size of model galaxies in each magnitude 
bin. Small gray dots are the observational data obtained by Courteau et al. 
(2007). 

5 Local galaxies 

In this section we compare the model predictions with local ob¬ 
servations. In what follows, we show the results of the u^GC- 
H2 model which has the highest mass resolution, unless other¬ 
wise mentioned. The adopted parameters related to the baryon 
physics are listed in table 2. 

5.1 Size and disk rotation veiocity of spirai gaiaxies 

First, we compare the predicted effective radius and disk rota¬ 
tion velocity of spiral galaxies with the observations. For the 
observational data, we use the data taken from Courteau et al. 
(2007). The sample of Courteau et al. (2007) is a compilation of 
the major samples of local spiral galaxies for which rotational 
velocities are available. Their sample Includes Mathewson et 
al. (1992), Dale et al. (1999), Courteau et al. (2000), Tully et 
al. (1996) and Verheijen (2001). The disk scale length of the 
sample galaxies are estimated from the /-band image, and the 
disk rotation velocities are estimated from Ha or H i line widths. 
Both of the disk size and the rotation velocity are corrected for 
inclination. 

Figures 7 and 8 show the scaling relations between the ef¬ 
fective disk radius and the /-band magnitude, and between the 
disk rotation velocity and the /-band magnitude (the so-called 
Tully-Fisher relation; Tully & Fisher 1977), respectively. The 
median and the 10th to the 90th percentile of the distribution 
of model galaxies in each magnitude bin are shown by black 
squares with error bars. The observational data are shown by 
small gray dots. The model very well reproduces these observed 
scaling relations over all magnitude range. The effect of the dy¬ 
namical response to the disk will be investigated in detail in 
future paper. 
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Fig. 3. The 1D (diagonal panels) and 2D (off-dlagonal panels) posterior probability distribution functions of five free parameters tuned in MCMC fitting. The 
solid vertical lines in diagonal panels show the median of each distribution. Probability distributions of all combinations of the five parameters are shown by 
grey contours. 
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Fig. 4. The r- and if-band LFs. The black solid line represents the model (i/^GC-SS) with the best-fit parameters determined by MCMC fitting. The black 
dotted line shows the model with the same parameters but using the high-resolution Ar-body data (i.'^GC-H2). The shaded region denotes the icr error in the 
model, estimated from the posterior probability distribution of the parameters (figure 3). The observational data shown in black filled circles are obtained by 
GAMA survey (Driver et al. 2012). 


Table 2. The list of free parameters related to astrophysical processes. Seven of these parameters, namely, r^^atar, Tatar, t4ot, 
“cool, esMBH> and Afaeedwere tuned to fit the local LFs and H i MF, by using a MCMC method. See text for details of parameter settings. 



Prior 

Posterior 

Median 10 to 90 percentile 

Meaning 

Q^star 

[-5.0, 0.0] 

-1.36 

[-1.14,-1.67] 

star formation-related 

S^star 

[0.01, 1.0] 

0.23 

[0.21, 0.26] 

star formation-related 


[0.0, 5.0] 

3.27 

[3.03, 3.52] 

SN feedback-related 

14ot (km/s) 

[50.0, 200.0] 

127.1 

[121.6, 133.1] 

SN feedback-related 

<3icool 

[0.1, 10.0] 

8.83 

[8.16, 9.55] 

AGN feedback-related 

IoSioC^smbh) 

[-2.0, 0.0] 

-0.50 

[-1.02, -0.14] 

AGN feedback-related 

loglo(A^seed/M©) 

[3.0, 6.0] 

5.45 

[4.83, 5.90] 

seed black hole mass 

/bh 


5 X 10“® (fix) 


fraction of the mass accreted onto SMBH 

during major merger 

rvo 

- 

2.5 X 10"® (fix) 

- 

coefficient of dust extinction 

/bulge 

- 

0.1 (fix) 

- 

major/minor merger criterion 

/mrg 

- 

1.0 (fix) 

- 

coefficient of dynamical friction time scale 

^diss 

- 

2.0 (fix) 

- 

energy loss fraction 


5.2 Size and velocity dispersion of elliptical galaxies 

In this subsection we compare the predicted effective radius and 
velocity dispersion of elliptical galaxies with the observations. 
For the observational data, we use the data compiled by Forbes 
et al. (2008). They take the central velocity dispersions of sam¬ 
ple galaxies from the catalog of Bender & Nieto (1990), Bender 
et al. (1992), Burstein et al. (1997), Faber et al. (1989), Trager 
et al. (2000), Moore et al. (2002), Matkovic & Guzman (2005) 
and Firth et al. (2007). The half-light radii are calculated from 
2MASS 7T-band 20th isophotal size, by using a empirical rela¬ 
tion based on Sersic light profiles (Forbes et al. 2008). 

Figures 9 and 10 show the scaling relations between the 
effective radius and the 7T-band magnitude, and between the 
velocity dispersion and the 7T-band magnitude (the so-called 
Faber-Jackson relation; Faber & Jackson 1976), respectively. 
The median and the 10th to the 90th percentile of the distri¬ 


bution of model galaxies in each magnitude bin are shown by 
black squares with error bars. The effective radius of model 
galaxy, re, is estimated from re = 0.744r6 (Nagashima & Yoshii 
2003), where ri, is the three-dimensional half-mass radius. The 
projected velocity dispersion is estimated as ctid = Vb/x/S after 
being increased to the central value by a factor of y/2 assuming 
the de Vaucouleurs profile. 

As shown in figures 9 and 10, our model underpredicts both 
of the size and the velocity dispersion of galaxies brighter than 
Mk ^ —20, comparing with the observations. The size and ve¬ 
locity dispersion are related to the dynamical mass of a galaxy 
as Mdyn oc TgVh, and therefore the model also underpredicts 
the dynamical mass of bright elliptical galaxies at a fixed mag¬ 
nitude. These results might imply that our treatment of bulge 
(and elliptical galaxies) formation process is oversimplified. We 
need to consider more realistic model for galaxy merger, as well 
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Fig. 6. The spatial distribution of the mock galaxies plotted on the past light-cone of an observer located at redshift zero. The color indicates the apparent 
magnitude of each galaxies in the 2MASS Ks-band. We only show the one-thousandth of galaxies which are randomly picked up from total sample, to avoid a 
confusion. 
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Fig. 8. /-band Tully-Fisher relation (i.e., disk rotation velocity against /-band 
magnitude) of spiral galaxies. The biack filled squares with error bars show 
the median and the 10th to the 90th percentile of the predicted disk velocity of 
model galaxies in each magnitude bin. Small gray dots are the observational 
data obtained by Courteau et al. (2007). 
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Fig. 9. Effective radius of elliptical and lenticular galaxies plotted against K- 
band magnitude. The black filled squares with error bars show the median 
and the 10th to the 90th percentile of the predicted size of model galaxies in 
each magnitude bin. Small gray dots are the observational data estimated 
by Forbes et al. (2008). 

as another channels of bulge formation such as disk instabili¬ 
ties. Furthermore, assumed IMF might also be responsible for 
the underprediction of mass-to-luminosity ratio. 

5.3 Cold gas 

Figure 11 presents the cold atomic hydrogen mass relative to 
the r-band luminosity against the r-band magnitude for local 
spiral galaxies. As described above, the atomic hydrogen mass 
of model galaxy is estimated as Mai = 0.54Mcoid (see section 
3.3). The median and the 10th to the 90th percentile of the 
distribution of model galaxies in each magnitude bin are shown 
by black squares with error bars. The observational data shown 
in small gray dots are taken from ALFALFA 40% catalog (a.40; 
Haynes et al. 2011). 

As already mentioned above, the uncertainties in the model 


Fig. 10. /T-band Faber-Jackson relation (I.e., projected central velocity dis¬ 
persion against /T-band magnitude) of elliptical and lenticular galaxies. The 
black filled squares with error bars show the median and the 10th to the 
90th percentile of the distribution of predicted velocity dispersion of model 
galaxies in each magnitude bin. The projected velocity dispersion of model 
galaxies are estimated as ctid = Vb/v^ after being increased to the cen¬ 
tral value by a factor of >/2 assuming the de Vaucouleurs profile. Small gray 
dots are the observational data compiled by Forbes et al. (2008). 

increase for galaxies having Hi mass less than 10®M© (see sec¬ 
tion 3.3). Furthermore, the a.40 catalog is highly incomplete 
for galaxies at Mhi < 10®M© (Haynes et al. 2011). Therefore 
we only plot galaxies having Mhi > 10®M© for both of the 
model and the observation. The diagonal solid line in figure 11 
corresponds to Mhi = 10®M©. 

We can see that the model very well reproduces the observa¬ 
tion over all magnitude range. The cold gas mass to luminosity 
ration is mainly determined by the balance of gas consumption 
rate by star formation and SN feedback. The agreement be¬ 
tween the model and observation seen in figure 11 supports the 
validity of our model of star formation and SN feedback. For 
more detailed discussion on the properties of cold gas compo¬ 
nent in our model, we refer the reader to Okoshi & Nagashima 
(2005) and Okoshi et al. (2010) although they are based on our 
previous model. 

5.4 Supermassive black holes 

In this subsection we present the properties of SMBHs at local 
universe. Figure 12 shows the predicted relation between the 
bulge mass and the SMBH mass, comparing with the observa¬ 
tional data obtained by McConnell & Ma (2013). To show the 
distribution of more massive and rarer objects, we also plot the 
iz^GC-M model in this figure. With a fixed mass fraction of cold 
gas accreted by a SMBH during a starburst (/bh = 0.005), the 
observed relation is naturally explained by the model. However 
/bh degenerates with other parameters which are related to 
bulge formation and SMBH evolution, such as /bulge and /mrg, 
and therefore we need another observational constraints to dis¬ 
cuss the physical meanings of these parameters. For example. 
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Fig. 11. Cold gas mass relative to r-band luminosity as a function of r-band 
magnitude for spiral galaxies. Smali gray dots are the observational data 
obtained by 40% catalog of ALFALFA survey (a.40; Haynes et al. 2011). 
Here we only show the galaxies having H i mass greater than io®Mq . The 
soiid diagonai line corresponds to the constant hydrogen mass of Mhi = 
10® Mq. The biack filled squares with error bars show the median and the 
10th to the 90th percentile of distributions of the model galaxies in each 
magnitude bin. We simply estimated the mass of cold atomic hydrogen as 
Mhi = 0.54Mcoid (see text for detail). 



Fig. 12. The SMBH mass versus bulge mass relation. The black filled 
squares with error bars show the median and the 10th to the 90 percentile of 
the distribution of model galaxies in each bin of bulge mass for the i/^GC-H2 
model. The crosses are the results of r/^GC-M model, which is shifted about 
+0.05 dex in log(Mbuige/M 0 ) to avoid a confusion. The observational 
data obtained by McConnell & Ma (2013) are shown by open circles. 

morphology-dependent LFs will help to resolve the degeneracy 
since /bulge and /mrg control the abundance of bulge compo¬ 
nent. Gravitational waves from SMBHs will also provide us 
strong and independent constraints (see, e.g., Enoki et al. 2004; 
Enoki & Nagashima 2007). Figure 13 show the predicted MF 
of local SMBHs, comparing with the observational estimation 
by Shankar et al. (2004). The model also well reproduces the 
observation over all mass range. 

For more detailed discussions on the properties of AGN 
populations, see Enoki et al. (2003), Enoki et al. (2014) and 
Shirakata et al. (2015), although they are based on our previous 
model. 
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Fig. 13. The mass function of local SMBHs. The analytical fit to the obser¬ 
vational data obtained by Shankar et al. (2004) is shown by the gray shaded 
region. The black solid line is the best-fit model. 

5.5 Distributions of gaiaxy coiors 

Figure 14 shows the distributions of (u — r) color of galaxies 
(i.e., differential number of galaxies per color bin) divided in 
several bins of the r-band magnitude. We compare the model 
predictions with the observed distributions extracted from Sloan 
Digital Sky Survey (SDSS) catalog (Baldry et al. 2004). As 
shown in figure 14, the model well reproduces the observed bi- 
modal distributions for the galaxies brighter than Mr = —19.5. 
However the model predicts systematically redder color for 
faint galaxies. This result might imply that the faint galaxies in 
our model obtain its stellar mass too early and have exhausted 
the almost all of cold gas, and consequently have redder colors. 
This discrepancy would be due to the oversimplified modeling 
of the star formation, SN feedback and stripping of hot gas in 
subhalos (cf. Makiya et al. 2014). We will investigate this issue 
in future paper. 

5.6 Main sequence of star-forming gaiaxies 

It has been known that the SFR and the stellar mass of star¬ 
forming galaxies are tightly correlated (the so-called “star¬ 
forming main sequence”; e.g., Brinchmann et al. 2004; Elbaz 
et al. 2007; Salim et al. 2007; Daddi et al. 2007). 

Figure 15 shows the SFR against the stellar mass for the 
model galaxies at z = 0. The star-forming galaxies and pas¬ 
sive galaxies are shown in blue and red dots, respectively. The 
black squares with error bars show the median and the 10th to 
the 90th percentile of the distributions of star-forming galaxies 
in each stellar mass bin. In the same figure we also show the 
observed relation obtained by Elbaz et al. (2007) by a solid line 
with typical errors by dashed lines. For the model galaxies, we 
adopt the same limiting magnitude, Mb < — 20 AB mag, with 
the sample of Elbaz et al. (2007). The separation criteria be¬ 
tween the star-forming galaxies and passive galaxies is also the 
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Fig. 14. Color distribution of galaxies (i.e., differential number of galaxies 
per color bin) in the different r-band magnitude bins (from top to bottom, 
-22.0 < Mr < -21.5, -20.0 < M, < -19.5, -18.0 < M, < -17.5 and 
-16.0 < Mr < -15.5). The biack solid lines in each panei are the analytical 
fit to the distribution of SDSS gaiaxies obtained by Baldry et al. (2004). The 
biack histograms are the model predictions. Both the model and observation 
are normalized to unity. 

same with Elbaz et al. (2007): the galaxies having (u— g) < 1.45 
are star-forming and the others are passive. Both of the SFR 
and the stellar mass are converted into those with Salpeter IMF 
from those with Chabrier IMF, by multiplying a factor of 1.8. 
We find that the model very well reproduces the observed tight 
correlation between the SFR and the stellar mass. 

5.7 Stellar-to-halo mass ratio 

Figure 16 presents the ratio of the stellar mass of central galaxy 
to the total baryon mass in their host halo against the total mass 
of their host halo. The total baryon mass Mbar is simply esti¬ 
mated as Mbar = Mh X (Ob/flm). TMs plot indicates an effi¬ 
ciency of star formation as a function of halo mass, and can be 
a tight constraint on the galaxy formation model. 

The median and the 10th to 90th percentiles of model galax- 



Fig. 15. The stellar mass vs SFR relation for local galaxies. Both the SFR 
and stellar mass are converted into those with Salpeter IMF from those with 
Chabrier IMF, by multiplying by a factor of 1.8. The solid and dashed lines are 
the observed relation and typical error obtained by Elbaz et al. (2007). The 
blue and red dots show the distribution of star-forming and passive galaxies 
in the model, respectively. Here we adopt the same color criteria with Elbaz 
et al. (2007), i.e., the galaxies having blue color ({u - g) < 1.45) are re¬ 
garded as the star-forming while the others are regarded as passive. For the 
model, we only plot the galaxies brighter than Mb = —20 AB mag, which is 
the limiting magnitude of the sample of Elbaz et al. (2007). The black filled 
squares with error bars show the median and the 10th to the 90th percentile 
of star-forming galaxies in each bin of stellar mass. 

ies in each halo mass bin are shown by the black squares with 
error bars. The solid and dashed lines show the average and 
Irr confidence level estimated by Moster et al. (2013) using an 
“abundance matching technique”, in which the halo mass is es¬ 
timated by matching the abundance of halos in Wbody simu¬ 
lations to the abundance of observed galaxies. The prediction 
of our model well agree with the result of Moster et al. (2013). 
The distribution of stellar-to-halo mass ratio has a peak around 
Mh ~ 10^^ Mq. It reflects effects of SN feedback and AGN 
feedback: the former efficiently works in lower mass halos be¬ 
cause the gravitational potential well is shallow in such halos, 
while the later efficiently works in massive halos because the 
cooling time is long enough and the central SMBH can suffi¬ 
ciently evolve in such halos. 

5.8 Mass metallicity relation 

Figure 17 shows the predicted relation between the stellar mass 
and the metallicity of cold gas for star forming galaxies. The 
median and the 16th to 84th percentile of the distribution of 
SDSS galaxies estimated by Tremonti et al. (2004) is also 
shown in figure 17 by solid lines. The cold gas metallic¬ 
ity is denoted by the gas-phase oxygen abundance in unit of 
12-Flog(0/H). The solar metallicity in this unit is 8.93 (Anders 
& Grevesse 1989). The metallicity with respect to the solar 
metallicity, Zq = 0.019 (Anders & Grevesse 1989), is also in¬ 
dicated on the right-side axis of figure 17 for reference. We de¬ 
fined the “star-forming galaxy” as a galaxy with specific SFR 
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Fig. 16. The stellar mass of central galaxies relative to the total baryon mass 
In their host halo as a function of the total mass of host halo. The black filled 
squares with error bars denote the median and the 10th to 90th percentile 
of model galaxies in each bin of the host halo mass. The black solid and 
dashed lines show the average and la confidence level of stellar mass ratio 
obtained by abundance matching technique (Moster et al. 2013). The total 
baryonic mass Mbar is estimated as Mbar = Afh x (ftb/STm). 



Fig. 17. Relation between stellar mass and cold gas metallicity, which is de¬ 
noted by the gas-phase oxygen abundance in unit of 12 -|-log(0/H). Solar 
metallicity in this unit is 8.93 (Anders & Grevesse 1989). The solid lines 
represents the 84th, 50th, and 16th percentile of local star-forming galax¬ 
ies observed by SDSS (Tremonti et al. 2004). The black filled squares with 
error bars show the 84th, 50th, and 16th percentile of the distributions of 
model galaxies in each magnitude bin. For the model, we defined star form¬ 
ing galaxy as a galaxy which are larger than 10“^^ yr~^ in specific star 
formation rate (i.e., SFR/M^tar) ■ 

(i.e., SFR/Mstar) higher than If we change this 

threshold to more lower value, for example, the relation will 
shift towards high-metallicity. 

Comparing with the observation, our model galaxies tend to 
have lower metallicities at the stellar mass range of Mstar < 
10^'^Mq. We will investigate this issue in future paper. 

6 Distant galaxies 

In this section we show the model predictions for the basic prop¬ 
erties of high -2 galaxies. 


6.1 Cosmic star formation history 

Figure 18 shows the redshift evolution of cosmic star forma¬ 
tion rate density (i.e., total SFR of all galaxies per unit co¬ 
moving volume). The blue solid line shows the result of our 
standard model (u^GC-Hl model). The SFR of model galaxies 
are converted into those with Salpeter IMF from those Chabrier 
IMF, by multiplying a factor of 1.8. The tx^GC-Hl model is 
also shown by red solid line to see the effect of mass resolu¬ 
tion. A discrepancy between these two models increases at high 
redshift, indicating that contributions from galaxies resides in 
lower mass halos become significant at high redshift. 

The standard model well reproduces the observations. At a 
redshift greater than 2 > 4, it seems that the model predictions 
are much greater than observed SFR density of Bouwens et al. 
(2014); however, their data only include galaxies brighter than 
M(1500A) < —17.0, while the other data and model predic¬ 
tions are integration in entire magnitude range. Furthermore, 
the survey of Bouwens et al. (2014) is designed to find galax¬ 
ies having blue colors, and therefore they might miss a popula¬ 
tion of dusty red galaxies. In fact, the predicted UV luminosity 
density (i.e., total luminosity of all galaxies per unit comoving 
volume) is roughly consistent with the data of Bouwens et al. 
(2014) when the effect of limiting magnitude are taken into ac¬ 
count (see next subsection). 

Our model predicts that a large amount of star formation ac¬ 
tivity has not yet been observed at distant universe. It will be 
investigated by future observations. 

6.2 Evolution of luminosity density in cosmic time 

Figure 19 shows the predicted redshift evolution of the lumi¬ 
nosity density at 1500 A(thick solid line). The intrinsic lumi¬ 
nosity density (i.e., without dust extinction) is shown by thick 
dotted line for reference. Note that the observational data plot¬ 
ted in figure 19 are not corrected for dust extinction effect 
thus they should be compared with the model with dust extinc¬ 
tion (thick solid line). As already mentioned above, the data 
of Bouwens et al. (2014) only include galaxies brighter than 
M(1500A) < —17.0. The model prediction taking into account 
the same magnitude limit with the Bouwens et al. (2014) is 
shown by thin solid line. We can see that the model well re¬ 
produces the observations. This result support a validity of our 
modelings of star formation and dust extinction. 

Figure 20 presents the redshift evolution of the sum of total 
IR luminosity (8-1000 ^m) of all galaxies per unit comoving 
volume. The total IR luminosity of model galaxies are esti¬ 
mated from the SED of each galaxy to be consistent with the 
total amount of stellar luminosity absorbed by dust. The obser¬ 
vational data are obtained by Gruppioni et al. (2013), by inte¬ 
grating the total IR LFs down to 10® L©. The model reproduces 
the observation within a factor of 2-3. The discrepancy between 
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Redshift z 

Fig. 18. The cosmic SFR density as a function of redshift. The red and 
blue solid iines represent the predictions by model with the W-body data of 
i/^GC-H1 (red) and !y^GC-H2 (blue), respectively. The parameters related 
to baryon physics are the same in these models. We also show the obser¬ 
vational data estimated by dust continuum (Pascale et al. 2009; Rodighiero 
et al. 2010; Karim et al. 2011) and UV continuum (Ouchi et al. 2004; 
Cucciati et al. 2012; Bouwens et al. 2014). The data of Hopkins (2004) are 
the compilation of various observations. All the data points are corrected for 
dust extinction, by the methods adopted in individual references. The data 
points of Bouwens et al. (2014) are obtained by integrating LF down to the 
Mab(1500 A) < -17.0, while the other observations and our model includes 
the contribution from all galaxies. The SFR of model galaxies are converted 
into those with Salpeter IMF from those Chabrier IMF, by multiplying a factor 
of 1.8. 



Redshift z 

Fig. 19. The redshift evolution of luminosity density at 1500 A. The filled cir¬ 
cles and open triangles are observational data compiled by Hopkins (2004) 
and obtained by Bouwens et al. (2014), respectively. The model prediction is 
shown by solid black line. For the purpose of comparison, we also show the 
model without dust extinction (dotted line). Those model predictions include 
a contribution from all galaxies. The data points of Bouwens et al. (2014) 
are obtained by integrating LF down to the Mab(1500A) < -17.0 while 
the other observational data are integrated down to zero luminosity. The thin 
solid line show the model prediction taking into account the magnitude limit 
of Mab(1500A) = -17.0. 


the model and observation is partly due to a contribution from 
AGNs, which is included in the observational data while not 
included in the model. 


6.3 Redshift evolution of A'-band luminosity function 

Figure 21 shows the redshift evolution of rest-frame A-band 
LF. The observational data are obtained by Cirasuolo et al. 
(2010). The model well reproduces the bright-end of LFs even 
at 2 = 2.0, which was not able to reproduce in our previous 
model. In new model, formation of massive galaxies are sup¬ 
pressed by AGN feedback only at low-redshift, and therefore 
the model can reproduce the bright-end LFs of local and high- 
2 galaxies at the same time. On the other hand, the model 
overestimates the abundance of dwarf galaxies over all redshift 
range. This discrepancy might suggest that SN feedback should 
be more efficient at high- 2 . However, there still remains some 
uncertainties in the observation. For example, cosmic variance, 
systematic error in fc-correction, and incompleteness of the sur¬ 
vey due to a surface brightness limit will affect the measurement 
of the faint-end slope of high -2 A-band LFs. 



Fig. 20. The model prediction for the redshift evolution of total IR luminosity 
density, comparing with the observational data Gruppioni et al. (2013). The 
total IR luminosity of model galaxies are calculated from the SED of each 
galaxy, to be consistent with the total amount of stellar luminosity absorbed 
by dust. 


7 Summary 

In this paper we present a new cosmological galaxy formation 
model, tz^GC, as an updated version of our previous model. 
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Fig. 21. The redshift evolution of rest-frame iT-band luminosity function. 
From top to bottom, we show the LFs at z = 0.5,1.0,1.5,2.0. The solid 
black lines are model predictions. The black filled circles with error bars are 
the observational data obtained by Cirasuolo et al. (2010). 

vGC (Nagashima et al. 2005; see also Nagashima & Yoshii 
2004). Major updates of the model are as follows: (1) the N- 
body simulations of the evolution of dark matter halos are up¬ 
dated (Ishiyama et al. 2015), (2) The formation and evolution 
process of SMBHs and the suppression of gas cooling due to 
the AGN activity (AGN feedback) is included, (3) heating of the 
intergalactic gas by the cosmic UV background is included, and 
(4) adopt a Markov chain Monte Carlo method in parameter tun¬ 
ing. Thanks to the updated A-body simulations, the minimum 
halo mass of the model reaches 1.37 x 10® Mq in the best case, 
which is below the effective Jeans mass at high redshift. In our 
largest simulation box (1.12 Gpc/h), we can perform statistical 
analysis for rare objects such as bright quasars. 

The main results of this paper are summarized as follows. 


1. We tuned the model to fit the local r- and JT-band LFs and 
HI MF by using a MCMC method. As a result, the model 
has succeeded well in reproducing these observables at the 
same time. 

2. The model well reproduces the scaling relations between 
the size and the magnitude, and the rotation velocity and 
the magnitude of spiral galaxies. For elliptical galax¬ 
ies, the model roughly well reproduces the observed size- 
magnitude relation and the velocity dispersion-magnitude 
relation. However, for bright elliptical galaxies, the model 
underpredicts both of the size and the velocity dispersion. 
We need to improve the model related to the galaxy merger 
and formation process of the bulge component. 

3. The model well reproduces the observed bimodal distribu¬ 
tion in color for bright galaxies. On the other hand, the 
model predicts redder color for dwarf galaxies comparing 
with the observations. This might be caused from our over¬ 
simplified prescription for star formation, SN feedback and 
stripping of hot gas. 

4. For massive galaxies (Mstar > 10^° M©), model well repro¬ 
duces the observed scaling relation between the stellar mass 
and gas phase metallicity at 2 = 0. However the model un¬ 
derpredicts metallicity of dwarf galaxies. This might also be 
caused by our oversimplified treatment of star formation and 
SN feedback. In addition, assumed IMF would also affect it. 

5. The observed scaling relation between the bulge mass and 
SMBH mass, and MF of local SMBHs are well reproduced 
in our model. 

6. The cosmological evolution of star formation rate density 
and UV luminosity density predicted by our model are well 
agree with the observations. We found that the model 
roughly reproduces the redshift evolution of total IR lumi¬ 
nosity density. We also compared the redshift evolution of 
the rest-frame JT-band LFs, and found that the model well 
reproduces the bright-end of LFs at 0 < 2 < 2. 

Since the main aim of this paper is to present the details of 
the calculation method of our model, we compared the model 
only with some basic observables mentioned above. Subsequent 
papers will discuss another topics related to galaxy formation: 
the clustering properties of quasars, the origin of cosmic NIR 
background, the properties of sub-millimeter galaxies, for ex¬ 
ample. 

The results of our model, including the LFs in several wave¬ 
bands, mass functions, and the mock galaxies, are publicly 
available on the web'. 
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